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ABSTRACT 

Propeller noise can be modeled as an amplitude modulated (AM) signal. 
Cyclic Spectral Analysis has been used successfully to detect the presence of 
analog and digitally modulated signals in communication systems. It can also identify 
the type of modulation. Programs for Signal Processing based on compiled 
languages such as FORTRAN or C are not user friendly, and MATLAB based 
programs have become the de facto language and tools for signal processing 
engineers worldwide. 

This thesis describes the implementation in MATLAB of two fast methods of 
computing the Spectral Correlation Density (SCD) Function estimate, the FFT 
Accumulation Method (FAM) and the Strip Spectral Correlation Algorithm (SSCA), to 
perform Cyclic Analysis. Both methods are based on the Fast Fourier Transform 
(FFT) algorithm. The results are presented and areas of possible enhancement for 
propeller noise detection and identification are discussed. 
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I. INTRODUCTION 


A. MOTIVATION 

Propeller related acoustic signatures typically exhibit modulation 
characteristics. These modulation characteristics originate from the cavitation 
process that takes place in the water due to the cyclic movement of the propeller. 

The cavitation process is basically the collapse of air and vapor bubbles 
due to variations in the static pressure. These variations in static pressure are a 
consequence of the passage of the propeller blades through the water. This 
movement, cyclic in nature, causes amplitude modulation in the static pressure, 
and as a consequence an amplitude-modulated (AM) signal can be detected in a 
receiver. 

Cyclostationary processing techniques have been used to detect and 
identify analog and digital communication signals very successfully. These 
techniques have the advantage of using a more realistic model for the signal 
than the stationary model used in most of the more conventional signal 
processing techniques. 

B. BACKGROUND 

The basic elements of cyclic spectral analysis are the time-variant cyclic 
periodogram and the time-variant cyclic correlogram. These two functions form a 
Fourier transform pair. This fact is known as the cyclic Wiener relation or the 
cyclic Wiener-Khinchin relation [Ref. 1: p. 49.] 
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In order to estimate the cyclic spectrum of a signal of interest, both the 
time-smoothed cyclic periodogram and the frequency-smoothed cyclic 
periodogram can be used, giving rise to two classes of computational algorithms: 
the time-smoothing algorithms and the frequency-smoothing algorithms. 

Although both classes of algorithms produce similar approximations to the 
cyclic spectrum, time-smoothing algorithms are considered to be more 
computationally efficient for general cyclic spectral analysis [Ref. 2: p.38]. Taking 
advantage of the efficiency, two computationally fast algorithms based on the 
time-smoothed cyclic periodogram were developed by Roberts et al. [Ref. 2]: the 
FFT Accumulation Method (FAM) and the Strip Spectral Correlation Algorithm 
(SSCA). 

C. THESIS GOALS 

The purpose of this thesis is to implement the FFT Accumulation Method 
as well as the Strip Spectral Correlation Method in MATLAB [ Ref. 3]. These two 
methods are already implemented in C [Ref. 4]. MATLAB is a more user-friendly 
and widely-used language that makes simulations and evaluations accessible for 
students and researchers. These cyclic spectral analysis programs written in 
MATLAB can easily be modified and are transportable across operating systems 
(i. e., UNIX, PC systems, MAC systems, VMS, etc). 

Both programs are used to compute the spectral correlation density (SCO) 
function estimate for a number of analog and digital modulated signals, such as 
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AM, PAM, ASK, and BPSK. The simulation results are then compared to the 
theoretical ones. After that the attention is focused on the AM signals for which a 
number of different signal-to-noise ratios (SNR) are processed using both 
methods. The results are presented and discussed for two types of modulation 
messages: a periodic message (a tone) and a random message (white noise). 

Results are presented by showing a three-dimensional plot of the cyclic 
spectrum surface, a contour plot, a plot of the power spectral density (PSD) 
obtained by setting the value of the cyclic frequency equal to zero, and some 
additional two-dimensional plots for cyclic frequencies of interest. 
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li. NOISE IN THE OCEAN 


Ross [Ref.5] defines noise as unwanted sound. The noise presence 
interferes with the signals that are of interest. If one is interested in detecting the 
presence of a particular class of surface ships, the sound generated by a near by 
submarine can be interpreted as noise. The reverse situation also leads to a 
similar statement. Therefore, the concept of noise has no absolute definition. It is 
a relative concept, and its characterization depends on the signal of interest in a 
particular situation. 

In this thesis, as we are interested in the detection and possible 
identification of the noise radiated by propeller, underwater noise is defined as 
any sound that interferes with our ability to detect and identify those signals. 

As reported by Ross [Ref 5: p. 1], undenvater noise is defined as sound 
in the water that limits the military effectiveness of naval detection and 
identification systems. Noise is unavoidable and sources that radiate as much as 
one watt of acoustics power can be detected at relatively long ranges by modern 
passive sonars. 

A. TYPES OF UNDERWATER NOISE 

There are several different sources of underwater noise that are grouped 
and classified in different ways according to the reference used. The main 
sources of underwater noise can be divided in ambient noise, self noise, and 
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radiated noise. Ambient noise and self-noise together constitute what is called 
the background noise which interferes with the operation of a sonar system. 

1. Ambient Noise 

Ambient noise is the prevailing, sustained unwanted background of sound 
at some location in the ocean. It excludes momentary, occasional sounds, such 
as the noise of a close-by passage of a ship or of an occasional rain squall. It is 
the background of noise, typical of the location and depth where a measuring 
hydrophone is located, against which a "signal," such as the sound of a 
submarine or the echo from a target, must be detected [Ref. 6: p. 1-1]. It 
comprises all noises associated with the medium in which a sonar operates that 
would exist in the medium if the sonar platform or vehicle itself were not present. 

The spectrum and characteristics of this kind of noise are complicated and 
depend upon location, position of the receiver, direction, and weather conditions. 
In its most general form, the ambient noise spectrum has some frequency bands 
where tonal components occur, and others where the spectrum is continuous 
and negatively sloping (“pink” noise), separated by portions where the spectrum 
is flat (“white” noise) or even reversed in slope [Ref. 6: p. 2-1]. This observation 
indicates that different sources of noise must exist and be prevalent in different 
regions of the spectrum. 

Urick [Ref. 6] divides the overall frequency range into five frequency 
bands: the ultra-low band {<\Hz), the infrasonic band (1 tolOHz), the low sonic 
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band {2Qto 200Hz), the high sonic band (200/o 50,000//z), and the ultrasonic 
{>50kHz). 

Almost nothing is known about the noise in the ultra-low band, since just a 
few measurements are reported. The infrasonic band contains the strong blade- 
rate fundamental frequency of propeller-driven vessels, plus one or two of its 
harmonics. It is of great interest to low frequency passive sonars. The low sonic 
band is characterized by the noise of distant shipping in areas where distant 
ships are prevalent. In areas remote from shipping lanes, the noise in this band 
is mainly dependent on wind speed [Ref. 6]. According to Ross [Ref. 5: p. 280], 
ship-generated noise is the dominant source of ambient noise in the spectral 
region between 20 and 200 Hz in most deep-water, open ocean areas and in 
highly traveled seas such as the Baltic and Mediterranean. The noise in the high 
sonic band is very dependent on the sea state and the wind speed. Thermal 
noise begins to dominate the noise background in the ultrasonic band. 

2. Self Noise 

Self noise is the noise associated with a platform and its sonar 
hydrophones. It includes the electronic noise generated by its preamplifiers, as 
observed by the sonar hydrophone array [Ref. 5: p. 4]. It can reach the receiver 
by transmission through the mechanical structure and by transmission through 
the water either directly from the source or by reflection from the sea surface. 

Self noise usually tends to increase as the speed of the platform 
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increases. At low frequencies and slow speeds, machinery noise dominates, 
whereas at high frequencies propeller and flow noise become important. 
Actually, as the speed is increased, these latter sources of noise assume more 
importance at all frequencies. 

At very low speeds, self noise is usually less important than ambient 
noise. At higher speeds the self noise can become the limiting factor [Ref. 7: p. 
413]. 

3. Radiated Noise 

Radiated noise is that sound radiated into the water which can be used by 
some receiver to detect the presence of the emitter at a considerable distance. 
It is very important for passive sonars, which are designed to exploit the 
peculiarities of this form of noise and to distinguish it from the background of self¬ 
noise and ambient noise in which it is normally observed [Ref. 8: p.328]. 

A detailed discussion of the mechanisms involved in the radiation of 
sound through the ocean can be found in several references such as Ross [Ref. 
5], Kinsler et al. [Ref. 7], and Urick [Ref. 8]. 

B. RADIATED NOISE FROM SHIPS, SUBMARINES, AND 
TORPEDOES 

The sources of noise on ships, submarines, and torpedoes can be 
grouped into three major classes: machinery noise, propeller noise, and 
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hydrodynamic noise. 

Machinery noise comprises that part of the total noise of the vessel 
caused by the ship's machinery. It originates as mechanical vibration of the 
many and diverse parts of a moving vessel. This vibration is coupled to the sea 
via the hull of the vessel. 

Machinery vibration can originate from five different sources. The first 
source of machinery noise are rotating unbalanced parts, such as out-of-round 
shafts or motor armatures. The second one are repetitive discontinuities, such as 
gear teeth, armature slots, and turbine blades. Reciprocating parts, such as the 
explosions in cylinders of reciprocating engines, are the third source of 
machinery vibration. The fourth are cavitation and turbulence in the fluid flow in 
pumps, pipes, valves, and condenser discharges. And mechanical friction, as in 
bearings and journals, is the fifth source of machinery noise. 

The first three sources of machinery vibration produce a line-component 
spectrum in which the noise is dominated by tonal components at the 
fundamental frequency and harmonics of the vibration-producing process; the 
other two give rise to noise having a continuous spectrum containing 
superimposed line components from structural members that are excited into 
resonant vibration. The machinery noise of a vessel may therefore be visualized 
as possessing a low-level continuous spectrum containing strong line 
components that originate in one or more of the repetitive vibration-producing 
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processes listed above. 

Even though the propeller is a part of the propulsion machinery of a 
vessel, the noise it generates has both a different origin and a different 
frequency spectrum from machinery noise. Machinery noise originates inside the 
vessel and reaches the \water by various processes of transmission and 
conduction through the hull. Propeller noise, on the other hand, originates 
outside the hull as a consequence of the propeller action and by virtue of the 
vessel's movement through the water. 

The source of propeller noise is principally the cavitation induced by the 
rotating propellers [Ref. 8: p. 334]. Because cavitation noise consists of a large 
number of random small bursts caused by bubble collapse, it has a continuous 
spectrum, covering a wide frequency range. 

Hydrodynamic noise originates in the irregular and fluctuating flow of fluid 
moving by the vessel. The noise created by the turbulent boundary layer is 
sometimes called "flow noise." Under normal circumstances, hydrodynamic noise 
is likely to be only a minor contributor to radiated noise, and is apt to be masked 
by machinery and propeller noises. 

C. PROPELLER NOISE 

Ross [Ref. 5: p. 253] describes cavitation of marine propeller as the most 
prevalent source of underwater sound in the oceans. Furthermore, when it 
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occurs, propeller cavitation is usually the dominant noise source for any single 
marine vehicle. Submarines and torpedoes often operate deep enough to avoid 
cavitation. Surface ships, on the other hand, generally have \well-developed 
propeller cavitation, with the result that their entire radiated spectra from as low 
as 5 Hz to as high as 100 kHz are controlled by this source [Ref. 5; p. 253]. 

Cavitation is essentially the rupture of bubbles in a liquid caused by 
reduction of local static pressure. It differs from boiling because the first is 
caused by a reduction of local static pressure while the second is due to an 
increase of temperature. Because of the pulse nature of the individual collapses 
and the random sequence of occurrence, the resultant spectrum covers a wide 
frequency range. Also, the pulsations of the aggregate volume of cavitation 
radiate strong tonals at frequencies below 70 Hz [Ref. 5: p. 270]. 

Of the various types of cavitation, blade-surface cavitation on the suction 
surface is the one that produces the highest noise levels. A more detailed 
explanation on the different types of cavitation, particularly on the blade-surface 
cavitation noise is found on Ross [Ref. 5: pp. 253-260]. 

Radiated spectra of surface ships are dominated by propeller cavitation 
noise except when the ships are operating at very slow speeds [Ref. 5: p. 272]. 
Some characteristics of surface ship noise that confirm the dominance of 
propeller cavitation are strong modulation of the broadband spectrum at shaft 
and blade frequencies and the radiation of low-frequency tonals at harmonics of 
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these frequencies. 

Propeller noise has been known for many years to be amplitude- 
modulated and to contain "propeller beats," or periodic increases of amplitude, 
occurring at the rotation speed of the propeller shaft, or at the propeller blade 
frequency, which is equal to the shaft frequency multiplied by the number of 
blades. The modulation at the shaft rotational frequency is due to slight physical 
differences between blades, that causes one blade to cavitate more than the 
others. It is this shaft-rate modulation that can be detected by the human ear and 
which enables an experienced sonar operator to determine the propeller 
rotational rate (rpm) and often classify the target [Ref. 5: p. 269]. Propeller beats 
have long been used by listening observers for target identification and for 
estimating target speed [Ref. 8: p. 338]. 

Propeller noise, with its origin in the flow of water about the propeller, 
creates tonal components in addition to the continuous spectrum of cavitation 
noise. Low-frequency spectra of cavitating ship propellers are usually dominated 
by tonal components at harmonics of the rotational frequency. More often, at the 
low-frequency end of the spectrum, propeller noise contains discrete spectral 
"blade-rate" components occurring at multiples of the rate at which any 
irregularity in the flow pattern into or about the propeller is intercepted by the 
propeller blades. The frequency of the blade-rate series of line components is 
given by 
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/„ = mns 

m 


(1) 


where /„ is the frequency, in Hertz, of the m'^ harmonic of the blade-rate series 
of lines, n is the number of blades on the propeller, and 5 is the propeller 
rotation speed in number of turns per second. 

Shaft and blade modulation frequencies for merchant ships are now 
significantly higher than they were during WWII. Forty years ago most merchant 
ships had three- or four-bladed propellers and operated at from 60 to 100 rpm. 
Shaft modulation frequencies were generally between 1.0 and 1.6 Hz and blade 
frequencies were from 3.5 to 6.5 Hz. Today, typical merchant propellers have 
four, five or six blades and operate at from 75 to 135 rpm; shaft frequencies 
range from 1.3 to 2.2 Hz, and blade frequencies are typically 6 to 12 Hz [Ref. 5; 
p. 279]. 
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III. CYCLOSTATIONARY PROCESSING 


The majority of the current signal processing techniques for intercepting or 
analyzing manmade signals in a noise contaminated environment typically 
utilize probabilistic models based on stationary statistics. In other words, they 
describe the signal on the average, and they have to restrict themselves to a 
small time interval in order for this approximation to hold. That limits the amount 
of data that can be used to derive the features in the signal and the resulting 
information. 

Most manmade signals, as are typically encountered in communication, 
radar, and sonar systems have some parameters that vary with time. Examples 
include sinusoidal carriers in amplitude, phase and frequency modulation 
systems, periodic keying of the amplitude, phase, or frequency in digital 
modulation systems, periodic scanning in some radar systems, and amplitude 
modulation in propeller noise. This requires that the random signal be modeled 
as cyclostationary, in which case the statistical parameters vary in time with 
single or multiple periodicity. 

Much of the background material in this chapter was taken from Gardner 
[Ref. 1]. 

A. CYCLOSTATIONARITY 

According to Gardner [Ref. 1], the essence of cyclostationarity is the fact 

that sinewaves can be generated from random data by applying certain 
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nonlinear transformations. As a consequence, a continuous signal x(/) is 


cyclostationary of order n (in the wide-sense) if and only if we can find some 
order nonlinear transformation of the signal, y{t)= f{x{t)), that will generate 
finite-amplitude additive sine-wave components, which produce spectral lines. In 
the same sense, a discrete-time signal x[m] is cyclostationary of order n (in the 

wide-sense) if and only if we can find some n'* order nonlinear transformation of 
the signal, y{rn\ = f{x[ni\], that will generate finite-amplitude additive sine-wave 
components, which will produce spectral lines [Ref. 1: p. 2]. 

A continuous signal y{t) contains a finite-amplitude additive sine-wave 
component with frequency a , a 0, if the Fourier coefficient 


M; ={yU)e-’^-") 


( 2 ) 


is not zero. In the same way, a discrete-time signal y[m] contains a finite- 
amplitude additive sine-wave component with frequency a ,a^0, if the Fourier 
coefficient 



( 3 ) 


is not zero. Here, denotes the sampling frequency and / is the square root of 
minus one. 

The operation (•) is the time-averaging operation defined as 
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( 4 ) 


(•) = lim — f (•) dt 


for analog signals, and as 


0 = 


1 


A'->C0 2N + 1 


Zo 


( 5 ) 


for discrete-time signals. 

For second-order cyclostationarity, the nonlinear transformation for a 
continuous signal x(r) (i. e., j(r) = /(x(r))) is given by 


y^{t) = x(t + T/2)x'(t-T/2), 

while for a discrete-time signal x[m] (i. e., }{m\ = /[x[w]]) it is given by 


( 6 ) 


yXm] = x[m]x‘[m-i\, (7) 

where the symbol * stands for complex conjugation. 

B. THE CYCLIC AUTOCORRELATION FUNCTION (ACF) 

The Fourier coefficient for second-order cyclostationarity is given by 

{x(t + r/2)x'(t-T/2)e‘'’^‘^ >. ( 8 ) 

This quantity is the fundamental parameter of second-order periodicity in 
continuous time and is called the cyclic autocorrelation function (ACF), R“{r), of 
x{t). 

For discrete-time signals, the ACF is defined as 
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(9) 


iTtai 

since delays other than sampling increments are not allowed. 

The ACF can be interpreted as measuring the amount of correlation 
between different frequency-shifted versions of a given signal, as shown in detail 
in Appendix A for an AM signal. 

Therefore, a signal exhibits second-order cyclostationarity in the wide- 
sense when its cyclic autocorrelation function, RU'^) for a continuous time 
signal or R“\l] for a discrete-time signal, is different from zero for some nonzero 
frequency a. The frequency a is called cycle frequency or cyclic frequency, and 
the discrete set of cycle frequencies a for which Rx(t)^0 or is called 

the cycle spectrum or cyclic spectrum. 

If a signal is cyclostationary, the cycle spectrum contains only harmonics 
(integer multiples) of the fundamental cycle frequency. Neverthless, if the signal 
has more than one fundamental cycle frequency, then the cycle spectrum 
contains harmonics of each of those frequencies. In this situation the signal is 
said to be polycyclostationary. 

C. THE SPECTRAL CORRELATION DENSITY FUNCTION (SCD) 

Signals usually have distinctive features in the frequency domain that are 
not easily seen in the time domain. Those features are generally used for 
detecting the presence of those signals. For instance, is very hard to detect the 
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presence of a sinusoidal signal, when embedded in noise, by just looking at its 
time-domain representation. The same signal can easily be detected in the 
frequency domain, provided the integration time can be made sufficiently long. 


For the same reason, it is beneficial to determine in the frequency domain 
the amount of correlation between frequency-shifted versions of x(r). The 
spectral correlation density function (SCD) is defined as the Fourier Transform of 
the cyclic autocorrelation function (ACF), to allow the localization in the 
frequency domain of the amount of time-correlation between frequency-shifted 
versions of the signal x(r). The SCD is given by 



( 10 ) 


for a continuous signal x(r), and by 


00 


«;(/)= Z KW'"’"' 


( 11 ) 


for a discrete-time signal x[«]. 

To determine the correlation in the frequency domain, we simply pass 
both of the two frequency translates «(/) = x(r)e'"'“' and v(r) = x(r)e'*“' through 
the same set of bandpass filters and then measure the temporal correlation of 
the filtered signals, as shown in Figure 1 (The signal «(/) represents a downshift 
in frequency by a/2 while v(/) represents an upshift in frequency by a/2 of the 
signal x(r)). In the limit when the bandwidth B of these narrowband components 
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approaches zero, we obtain 


5r(/)=limi([/.;(()0t-O)](^l(')®v(O]' ). (12) 

where the symbol 0 stands for convolution, and h{{t) is the impulse response 
of the bandpass filters. 



Figure 1 Spectral Correlation Analyzer (after Gardner [Ref. 1]) 


Strictly speaking, the SCO is not a valid density function in the usual 
sense, since it is not nonnegative and, in fact, not even real-valued. However, 
because of the similar properties that the SCO does share with the PSD, the 
term density is retained. 
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The SCD is also called the cyclic spectral density. Observe that although 
the argument / of the SCD is continuous, as it always will be for a random 
signal, the argument a is discrete, as it always will be since it represents the 
harmonic frequencies of periodicities underlying the random time-series. 

A detailed example of the computation of the SCD for an amplitude- 
modulated (AM) signal is given in Appendix A. 
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IV. ESTIMATION OF THE SPECTRAL CORRELATION DENSITY 

FUNCTION 


Cyclic spectral analysis is used to detect the presence of a signal via the 
spectral correlation density (SCD) function. To accomplish this goal a series of 
codes that estimates the SCD function were developed in MATLAB language. 
Those codes are implementations of two FFT based time smoothing algorithms 
called the FFT Accumulation Method (FAM) and the Strip Spectral Correlation 
Algorithm (SSCA). The majority of the background material in this chapter was 
taken from Roberts et al. [Ref. 2]. 

As reported by Roberts et al. [Ref. 2], cyclic spectral analysis algorithms 
generally fall into two classes: those that average in frequency (frequency 
smoothing) and those that average in time (time smoothing). Although both 
classes of algorithms produce similar approximations to the cyclic spectrum, time 
smoothing algorithms are considered to be more computationally efficient for 
general cyclic spectral analysis. 

Both methods are based on modifications of the time smoothed cyclic 
cross periodogram, defined as: 




(13) 


where Ar represents the data time span, and X-An,/ + —I and Tj 


a 


nJ-\ 


are 


the complex envelopes of narrow-band, bandpass components of the signals 
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x(n) and y{n), respectively. The complex envelopes are also called the complex 
demodulates of the signals. These complex demodulates are computed in the 
following way: 

N-/2 

(14) 

r=-N'l2 

and 

N'f2 

yA«J)= (15) 

r=-N'f2 

where a(r)is a data tapering window of length T = N'T^ seconds, with being 
the sampling period. The Fourier transform of plays the role of a spectral 
window. The particular shape of window, especially the spectral window, is of 
considerable importance. Windows other than the rectangle have a tapering 
effect on the data they multiply, since data occurring away from the aperture 
center are attenuated relative to the data at the aperture center. A data tapering 
window whose Fourier transform has low skirts and low sidelobes is desirable to 
reduce cycle leakage. A Hamming window is used for the input bandpass filters. 

Figure 2 shows a basic implementation of the discrete time smoothed 
cyclic cross periodogram, where the symbol * stands for complex conjugation. 
The complex demodulate frequencies /, and are related to the spectrum 
frequency / and the cyclic frequency a of the point estimate by the following 
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two equations; 


/l +/2 

f = (16) 

and 

cc = f,-h. (17) 



Figure 2 Discrete Time Smoothed Cyclic Cross Periodogram (after Roberts et al. 
[Ref. 2]). 


A. FFT ACCUMULATION METHOD (FAM) 

In this method, the complex demodulates are estimated efficiently by 
means of a sliding -point FFT, followed by a downshift in frequency to 
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baseband. In order to allow for an even more efficient estimation, the //'-point 
FFT is hopped over the data in blocks of L samples (channelization). This 
means that L data points are skipped between successive computations of the 
//'-point FFT. The value of L was chosen to be equal to //'/4 since, according 
to Reference 2, p.462, it allows for a good compromise between maintaining 
computational efficiency and minimizing cycle leakage and cycle aliasing. The 
value of //' is determined according to the desired resolution in frequency (A/) 
used in the algorithm, and is given by 

N'=^. (18) 


//' is chosen to be the power of 2 equal to or larger than the number 
given by Eq. 18 to take advantage of the Fast Fourier Transform (FFT) algorithm 
without making use of zero-padding. 

After the complex demodulates are computed and the product sequences 
between each one of them and the complex conjugate of the others are formed, 
the time smoothing is accomplished by means of a P-point FFT. The value of P 
is determined according to the desired resolution in cyclic frequency (Aa), and is 
given by 


P = 


fs 

L Aa 


(19) 


Again, P is chosen to be the power of 2 equal to or larger than the 
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number given by Eq. 19 to take advantage of the FFT algorithm avoiding zero¬ 
padding. 

Figure 3 illustrates the generation of the complex demodulates while 
Figure 4 shows the implementation of the FAM method. 





Xr(kL,f) 


Figure 3 Computation of the Complex Demodulates 





Xj 






P-point FFT 



Figure 4 Implementation of the FFT Accumulation Method 
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An advantage of having complex demodulates is that there is no need to 
worry with a correction factor to take care of the phase shift introduced by 
overlap processing. The last multiplier in Figure 3 (i. e., complex exponential) 
provides the correction to compensate for the overlap processing artifacts. 

The MATLAB codes that compute and plot the SCO function estimate 
using the FAM method are called AUTOFAM and CROSSFAM. AUTOFAM 
computes and plots the auto spectral correlation density function (amount of 
correlation between frequency shifted versions of a given real or complex valued 
signal) estimate. CROSSFAM computes and plots the cross spectral correlation 
density function estimate for two different real or complex valued signals. 

The inputs required for these programs are the signal(s), the sampling 
frequency (/J, the desired frequency resolution (A/), and the desired resolution 
in cyclic frequency (Aa). In the case of two signals, the sampling frequencies 
must be the same. 

The programs are listed in Appendices B and C. 

B. STRIP SPECTRAL CORRELATION ALGORITHM (SSCA) 

In the SSCA algorithm, the complex demodulate of one of the signals is 
computed in the same way it is done for the FAM method (channelization). The 
complex demodulated sequence is directly multiplied by the complex conjugate 
of the other signal. Then, the resultant signal is smoothed in time by means of an 
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N -point FFT. Here, N is the total number of data points {N = PL). 

Figure 5 shows the implementation of the SSCA. The complex 

demodulated component Xj[kL,f^^ is obtained as shown in Figure 3. 

It appears that the complex demodulate ) is at a sampling rate 

fJL due to the introduction of the decimation factor L , and consequently it can 
not be directly multiplied by which is at a sampling rate . However, the 
demodulated sequence is interpolated to match the sampling rate of x[«] by 
means of a process called replication. Replication is accomplished by holding the 
value of each complex demodulate sample for L samples [Ref.2]. 

The MATLAB codes generated to compute and plot the SCO function 
estimate for given signal(s) using the SSCA method are called AUTOSSCA and 
CROSSSSCA. AUTOSSCA computes and plots the auto spectral correlation 
density function (amount of correlation between frequency shifted versions of a 
given signal) estimate. CROSSSSCA computes and plots the cross spectral 
correlation density function estimate for two real or complex valued signals. 

The inputs required for these routines are the signal(s), the sampling 
frequency (/J, the desired frequency resolution (A/), and the desired resolution 
in cyclic frequency (Aa). In the case of two signals, the sampling frequencies 
must be the same. 

The codes are given in Appendices D and E. 
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Figure 5 Implementation of the SSCA 









V. EXPERIMENTAL RESULTS 


In this chapter, the theoretical results of computing the SCD function for 
some analog and digital modulated signals are presented, together with the 
experimental results obtained as output from the programs AUTOFAM and 
AUTOSSCA. An interpretation of the results is also provided. 

A. ANALOG MODULATED SIGNALS 

1. Amplitude-Modulated (AM) Signal 

Consider the following amplitude-modulated (AM) signal, x[«], given by 

x{n] = a{n]f{n \, (20) 

where a[n] is a purely stationary low pass message signal with power spectral 
density 5'^(/), and /?[«] is a sinusoidal carrier wave given by 

p\n] = cos(2;r/„« + (21) 


In Eq. 21, /^ = F^jF^ is the carrier digital frequency and <j)^ is the carrier 
phase. F^ is the carrier frequency in Hz and F^ is the sampling frequency in 
Hz. 

The SCD function for this amplitude-modulated (AM) sine wave is given 
by 


31 


sf(/)= 


^S,{/+/J +jS,(/-/J. « = o 

I s, (/)£“«-, a = ±2/. 

0, otherwise. 


( 22 ) 


Details of the derivation are given in Appendix A. 

When a[n] is a tone, some of the assumptions made to obtain Eq. 22 are 
not valid. If the message a[n\ is a tone at digital frequency /„ given by 

a[«] = cos(2;r/,«), (23) 

then it is necessary to go through the same steps as for a stationary lowpass 
signal with no spectral lines in the message’s power spectral density (PSD), as 
shown in Appendix A. 

Let us assume that x[«] is an amplitude-modulated single-sideband 
(AMSSB) signal, obtained by transmitting only the lower-side frequency, given by 

x[«] = ^[cos(2;r/„«) • cos(2;r/,« + <f>„) + sin(2;r/„«) • + <f >^)] 

= ^ cos[2;r(/„ -/„)« + ^„ ]. (24) 

Since the cyclic autocorrelation function (ACF) is given by 

RM = (4«] An - . (25) 

replacing x[n] and x[n-i] by Eq. 24 leads to. 
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j_^w[a-2(/„-/,)]«^^/2;r[o-2(/„-/.)]«1 ^^o+2(/.-/„)]f^^-/2;r[a+2(/<,-/„)]«^^-,2?S„ ^20) 

This can be rewritten as 


R:[e]= 


^cos[2;r(/„-/j£], a = 0 

« = ±2(/„-/„). 




(27) 


The SCD is the Fourier transform of the cyclic autocorrelation function. 
Therefore, the SCD is given by the Fourier transform of Eq. 27 leading to 




^Kz-z.+zJ+^z+z.-/,)], «=o 


7^(Z>*'^'-. 


a = ±2(/.-/.). 


(28) 


Now, let us suppose that x[n] is an amplitude-modulated double¬ 
sideband supressed carrier (AMDSB-SC) signal given by 

x[n] = cos(2;r ■ cos(2;r f^n + ^^) • (29) 

This can be written using a trigonometric identity as 

^ {cos[2;r(X, + ] + cos[2;z-(/„ -/„)« + ^„ ]} • (30) 

So, replacing x[«] and x[n-i\ by Eq. 30 into Eq. 25 leads to. 


^;W = |[cos2;r(/, +/„)f+ cos2;r(/, -/,)f]e“(e-"'“> + 


J_ ^i>r[a-2(/„+/J]« /-(2s-[a-2(/„+/,)]» \ J_ ^/;r[a+2(/„-/J]f l-i2iz{a-2f^)n \ 

16 ' / 16 ' ' 




16 
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l_^;*[o+2(/,-/„)]f^^-;2;r(a+2/„)n\^^-,2«l„ _L ^'^[«-2(/„-/<,)]^/^-<2;r(a+2/Jn 


16 


e H- e 

16 


^^-<2;r(a+2/Jny 


]_^/;r[o+2(/„+/.,)]<^^-/2;r(a+2/Jn^ + _LgM“-2(/»+/»l]^^g-'2'r(a-2/„)«^^/2<it„ 


16 


16 


J_^/a-(a-2(/<,-/<,)]f/^-*2;r[a-2(/,-/<,)]«\^/2<H„ ^ J_^w[a+2(/<,+/Jlf ^g-/2;r(ar+2/<,)n ^ 


16 


16 


J_ ^w[a-2(/„+/„)]? ^^-/2;r(a-2/,)nj_^J_^/;r[a+2(/„-/<,)]? ^^-;2;r|a+2(/„-/<,)]« ^^-> 2 «(„ 


16 


16 


This can be rewritten as 


R:[e]= 


^[cos2;r(/, + f^)e + cos2;z-(/„ - f^)l\ a = 0 

« = ±2/<, 


— j. -1 


g—v<,* + —e =-cos2;^ 


16 16 


,±'2«><. 


1 


16 


coslTifJ, 


j_ 

16' 


_g±'2«><. 


a=±2(/.-/,) 

«=±2/. 
a = ±2(/.+/,). 


(32) 


Therefore, the SCD is given by the Fourier transform of Eq. 32 leading to 


5;(/)= 


:|^K/ - /. - L) + <?(/-/<. +/.) fo )+^/ +/. + /.)], 

a=0 


« = ±2/. 


a = ±2(/„-/„) 


« = ±2/. 


a = ±2(/„+/J. 


(33) 


For the case where x[a 2 ] is an amplitude-modulated double-sideband 
transmitted carrier (AMDSB-TC) signal given by 

x[«] = [l + cos(2;r/„«)] • cos(2;r/„« + (j)^). (34) 
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This can be written using a trigonometric identity as 


x[n] = cos(2;7’ f^n ) + ^{cos[2;r(4+/J « + ^^] + CO s[2;r(/^+ (35) 


So, replacing x[n\ and x[n-^] by Eq. 35 into Eq. 25 leads to, 

K [<] = {|cos(2;r//) + |[cos 24 /, +/,)£ + cos2;r(/, - /„)4 ) + 




+ -e 
8 


+ -e ' +-e 

8 8 


+ 


J_ ia-(a+2/,)l/ -;2a-(a-/„)n 

l8 8 8 8 |\ 


I J_ w[a- 2 (/,-/„)]^ , _L w[o+ 2 (/.+/.)]^ 1 / -/ 24 «+ 2 /Jn 
1 16" ■'16" )\" 


+ 




J_ w[or+ 2 (/,-/„)]^_/ 2 «l/-; 24 a+ 2 (/„-/„)]» 

16 \ 


e ’"“(e 


16 




lgMo+2/„)f |^-;2«l„ ^^-/ 24 a:+( 2 /„-/J]n 


l^lg'^(«+2/„)e ^ J_ M«+2(/„-/J]f J_^;>[o+2(/,+/.)]^l /^-»2^a+2/Jn\ 

[4 16 16 J ^ ^ 

|lgM“-2/.)^ + J_g^“-2U-/o))^ ^ _LgM“-2(/»+/«)l^lg'2«(„ r-/2ff(a-2/„)«\ 

[4 16 16 J ' ^ 


35 


j_ ,>|a+2(/„+/o)]f 1 /;r(o+2/„)H -/2(0„ / -' 2 ;r[a+( 2 /„+/„)]« 


+ — e 
8 


+ 


4. L''24 /„-''2'f[«-(2/„+/„)]« 


+ 


J_ ^ i)r\ a+2{f^ +/«)] (^ -/2;r|a+2(/„ +/„)]« 


16 


J_ w[a-2(/„+/„)]f^/2j>„ / -i2;r[a-2(/„+/<,)]« 

16 \ 


( 36 ) 


This can be rewritten as 




|Cos[ 22 r(/„ -f^y\ + ^coi27tfj) + ^co^n{f^+ f^)i\ a = 0 

j cos[;r(2/„ - /. )^] + ^ cos[ 2 r{ 2 /„ + /, )4 “ = -fa 

a = ±2f 


\cQ 42 nfj), 


J_ 

16 ' 




-cos(2r/j)e*'''^% 
[7 + ^cos(2;r/,£) 


±/ 2 <^o 


1 

-< 

16 




« = ±2(/.-A)(37) 
« = ±(2/„-/J 
« = ± 2 /„ 

« = ±(2/„+/J 
« = ±2(/„+/J. 


Therefore, the SCD is given by 

-±/. ( 38 ) 

a = i2f, 
a = n[f,-f,) 


s;(/)= 
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“ = ±( 2 /,-/,) 
a = ±2/„ 

a = ±(2X-/,) 

a = ±2(/, + /,). 
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The results obtained from AUTOFAM and AUTOSSCA are presented as a 
sequence of three plots, for each of the signals utilized. This presentation is 
maintained throughout (i. e., Figure 6 - Figure 41). The first plot is a surface plot 
that shows the magnitude of the SCD function estimate in each region of the 
bifrequency plane with coordinates / and a. The second plot is a two- 
dimensional contour plot. It gives a better view of the position of the features in 
the bifrequency plane. The third plot is a set of two-dimensional slices of the 
SCD function estimate for fixed values of the cyclic frequency a . 

Figures 6-11 show the results obtained by using the programs 
AUTOFAM and AUTOSSCA on an amplitude-modulated single-side band 
(AMSSB) signal. For f^=5\2Hz and /,=2048//z, Equation 28 leads to the 
following result: 



a = 0 


(39) 


a = ±3072ife. 


So, according to Eq. 39, ones expect to obtain peaks at / = ±\536Hz for 
a = 0, and at f = 0 for a = ±3012Hz. The results in Figures 6-11 are in 
agreement with the theoretical results in Eq. 39 and show that AUTOFAM gives 
a clearer picture than AUTOSSCA. This is even true when AUTOSSCA uses a 
finer resolution in cyclic frequency (i. e., Aa = 64Hz for AUTOFAM versus 
Aa = 32 Hz for AUTOSSCA). 

The results obtained for an amplitude-modulated double-side band 
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suppressed-carrier (AMDSB-SC) are presented in Figures 12-17. For 
/„ = 5\2 Hz and = 2048//z, Equation 33 leads to the following result: 


■ 5 . 1 /) = 


^[<5{/ - 2560)+<5(/ -1536)+<5(/ +1536)+ S{f +2560)], a = 0 

^[<5{/-2048)+^/+2048)], a = ±1024/fe 

—^/), a = mi2Hz 

16 

^[<5{/-512)+<5{/+512)], a = ±m6Hz 

’ ' ' a = ±Sl20Hz. 




(40) 


According to Eq. 40, ones expect to have four peaks at / = ±l536Hz and 
/ = ±25607fe, for a = 0 ; two peaks at / = ±2048//z, for or = ±1024//z; one peak 
at / = 0, for a = ±3072Hz] two peaks at / = ±5\2Hz, for a = ±4096Hz‘, and a 
peak at / = 0 , for a = ±5\20Hz . Figures 12-17 confirm the theoretical results as 
given by Eq. 40. 

The results for an amplitude-modulated double-side band transmitted- 
carrier (AMDSB-TC) are presented in Figures 18-23. For f^=5\2Hz and 
= 294% Hz , Equation 38 leads to the following result; 


5 . 1 /)= 


■i.5{/+2560)+ii;(/+2048)+/^(/ + 1536)+.lrf{/-I536) + lrf(/-2048)+/<5{/-2560). a = 0 

l<5(/ + 2304)+l^/ + 1792)+i<5(/-1792)+^^/-2304), a = ±5127* 

■/<S(/+2048)+/<5(/-2048). a = ±10247/: 

a = ±307277: 

\-e{f + 256)+i <!>(/ - 256) , a = ±358477: 

|_8 8 

j^.i.5(/+5l2)+l7i(/)+-l<5{/512.) a = ±409677: 

\-S[f+ 256)+- ^/ - 256) , a = ±4608//r 

[_8 8 

a = ±512077:. 

16^ ^ 


(41) 


According to Eq. 41, ones expect to have two big peaks at / = ±2048//z 
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and four smaller peaks at / = ±1536//r and / = ±2560//r, for a = 0; four peaks 
at / = ±\192Hz and / = ±2304//z , for a = ±5\2Hz \ two peaks at / = ±2048//z, 
for O' = ±1024/yz; one peak at / = 0, for o' = ±3072//z; two peaks at / = ±256Hz , 
for O' = ±35MHz ; a big peak at / = 0 and two smaller peaks at f = ±5\2Hz, for 
a = ±4096Hz ; two peaks at f = ±256//z, for o' = ±4608//z; and a peak at / = 0, 
for a = ±5\29Hz. Figures 18-23 verify the theoretical results of Eq. 41. 

To have a reliable estimation of the spectral correlation density function is 
necessary that the product A/A/ »1 [Ref. 1]. This condition requires that 
A/ »Ao'. In some of the results obtained (i. e.. as in Figure 21). this 
requirement was not met. Computational limitations did not allow for a better 
resolution in the plots, since this translates to generating a large amount of data, 
making impossible to obtain a plot on the printer and/or on the screen with the 
available PC and workstation resources. 
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Figure 6 Surface plot of the SCD estimate magnitude for an Amplitude 
Modulated (AM) Single-Side Band Signal, using AUTOFAM, with the following 
parameters; Af = 256Hz, Aa = 6AHz, f^ = 204ZHz, /, =512i/z, and /^ =8192/fe, 
where is the carrier frequency, /, is the tonal frequency, and is the 
sampling frequency. 
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Figure 7 Contour plot of the SCO estimate magnitude for an Amplitude 
Modulated (AM) Single-Side Band Signal, using AUTOFAM, with the following 
parameters: A/ = 256/*, ha = 64/*, = 2048/*, /, =512/*, and /, = 8192/*, 

where is the carrier frequency, /, is the tonal frequency, and is the 
sampling frequency. 
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Figure 8 Plots of the SCD estimate magnitude for an Amplitude Modulated (AM) 
Single-Side Band Signal, using AUTOFAM, for a = 0 and a = 3072Hz, 
respectively, vwith the follovving parameters; A/ = 256Hz, Aa = 64Hz, 
/. = 2048/fe, /, =5\2Hz , and = S192Hz, where is the carrier frequency, /, 
is the tonal frequency, and is the sampling frequency. 
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Figure 9 Surface plot of the SCD estimate magnitude for an Amplitude 
Modulated (AM) Single-Side Band Signal, using AUTOSSCA, with the following 
parameters: A/ = 256/fe, Aa = 32/fe, /, =2048//z, f,=5\2Hz, and /, = 8192//z, 
where /. is the carrier frequency, /, is the tonal frequency, and 's the 
sampling frequency. 
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Figure 10 Contour plot of the SCO estimate magnitude for an Amplitude 
Modulated (AM) Single-Side Band Signal, using AUTOSSCA, with the following 
parameters: A/ = 256Hz, A« = 32//z, /^ = 2048//z, /, =512//z, and /^ = 8192//z, 
where is the carrier frequency, /, is the tonal frequency, and f, is the 
sampling frequency 


44 



























Figure 11 Plots of the SCD estimate magnitude for an Amplitude Modulated 
(AM) Single-Side Band Transmitted Carrier Signal, using AUTOSSCA, for a = 0 
and a = 3072Hz , respectively, with the following parameters; A/ = 256Hz, 
Aa = 32Hz, /, = 2048/fe, f,=5\2Hz, and f^ = Sl92Hz, where /, is the carrier 
frequency, /, is the tonal frequency, and /,'s the sampling frequency. 
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Figure 12 Surface plot of the SCD estimate magnitude for an Amplitude 
Modulated (AM) Double-Side Band Supressed Carrier Signal, using AUTOFAM, 
with the following parameters; A/ = 256Hz, Aa = 7>2Hz , = 2048//z , f, =5\2Hz , 
and /, = 8192//z, where is the carrier frequency, /, is the tonal frequency, 
and is the sampling frequency. 
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Figure 13 Contour plot of the SCD estimate magnitude for an Amplitude 
Modulated (AM) Double-Sided Band Supressed Carrier Signal, using 
AUTOFAM, with the following parameters: A/ = 256Hz , Aa = 32i/z, /<. = 2048//z, 
/, =5\2Hz, and /^ = 8192//z, where is the carrier frequency, /, is the tonal 
frequency, and is the sampling frequency. 
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Figure 14 Plots of the SCD estimate magnitude for an Amplitude Modulated 
(AM) Double-Side Band Supressed Carrier Signal, using AUTOFAM, for a = 0, 
a = \024Hz, a = 3012Hz, a = 4096Hz, and a = 5\20Hz, respectively, with the 
following parameters: Af = 256Hz, Aa = 32Hz, /. = 2048//z, f,=5\2Hz, and 
= 8192//z, where /, is the carrier frequency, /, is the tonal frequency, and 
is the sampling frequency. 
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Figure 15 Surface plot of the SCD estimate magnitude for an Amplitude 
Modulated (AM) Double-Side Band Supressed Carrier Signal, using 
AUTOSSCA, with the following parameters: Af = 256Hz, Aa = 22Hz, 
= 2048/fe, /, =S\2Hz , and = %\92Hz , where is the carrier frequency, /, 
is the tonal frequency, and is the sampling frequency. 
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Figure 16 Contour plot of the SCO estimate magnitude for an Amplitude 
Modulated (AM) Double-Sided Band Supressed Carrier Signal, using 
AUTOSSCA, with the following parameters: A/ = 256/fe, Aa = 32Hz, 
/.=2048//z, /, =512/fe,and /, =8192/fe, where is the carrier frequency, f, 
is the tonal frequency, and is the sampling frequency. 
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Figure 17 Plots of the SCD estimate magnitude for an Amplitude Modulated 
(AM) Double-Side Band Supressed Carrier Signal, using AUTOSSCA, for a = 0, 
a = l024Hz, a = 3072Hz, a = 4096Hz, and a = 5120i*, respectively, with the 
following parameters: A/ = 256Hz, Aa = 32Hz, = 204%Hz, f, = 5\2Hz, and 

f, = Sl92Hz, where is the carrier frequency, /, is the tonal frequency, and 
is the sampling frequency. 


51 




































alpha(Hz) 


Figure 18 Surface plot of the SCD estimate magnitude for an Amplitude 
Modulated (AM) Double-Side Band Transmitted Carrier Signal, using AUTOFAM, 
with the following parameters: A/ = 128//z, Aa = 32Hz, f^ = 204SHz, f,=5\2Hz, 
and f^ = S\92Hz, where is the carrier frequency, /, is the tonal frequency, 
and is the sampling frequency. 
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Figure 19 Contour plot of the SCO estimate magnitude for an Amplitude 
Modulated (AM) Double-Sided Band Transmitted Carrier Signal, using 
AUTOFAM, with the following parameters: A/ = 128/fe, Aa = 32//z, = 2048//z, 

/, =512^fe, and = %\92Hz, where is the carrier frequency, /, is the tonal 
frequency, and is the sampling frequency. 
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Figure 20 Plots of the SCD estimate magnitude for an Amplitude Modulated 
(AM) Double-Side Band Transmitted Carrier Signal, using AUTOFAM, for a = 0, 
a = 1024/fe, a = 2(yj2Hz, a = 4096Hz, and a = 5\20Hz, respectively, with the 
following parameters; A/ = 128/fe, Aa = 32Hz, f^=2048Hz, f,=5\2Hz, and 
= 8192/fe, where is the carrier frequency, /, is the tonal frequency, and 
is the sampling frequency. 
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Figure 21 Surface plot of the SCD estimate magnitude for an Amplitude 
Modulated (AM) Double-Side Band Transmitted Carrier Signal, using 
AUTOSSCA, with the following parameters: A/ = 128i/z, Aa = 128//z, 
= 2048//z, /, = 512/fz, and /, = 8192//z , where is the carrier frequency, /, 
is the tonal frequency, and is the sampling frequency. 
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Figure 22 Contour plot of the SCO estimate magnitude for an Amplitude 
Modulated (AM) Double-Sided Band Transmitted Carrier Signal, using 
AUTOSSCA, with the following parameters; A/ = 128//z, Aa = 32//z, 
X = 2048/fe, /, =5\2Hz , and = 8192//z, where is the carrier frequency, f, 
is the tonal frequency, and is the sampling frequency. . 
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Figure 23 Plots of the SCD estimate magnitude for an Amplitude Modulated 
(AM) Double-Side Band Transmitted Carrier Signal, using AUTOSSCA, for 
a = 0, a = 1024/fz, a = 3012Hz, a = 4096Hz, and a = 5120/fe , respectively, with 
the following parameters; A/ = 32Hz, Aa = 32Hz, f^ = 204%Hz, f,=5\2Hz, and 
= 8192//Z , where is the carrier frequency, /, is the tonal frequency, and f, 
is the sampling frequency. 
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2. Pulse-Amplitude Modulated (PAM) Signal 


In the previous section, consider the following form for the carrier wave, 

p[n], 

pM = X ) (42) 

OT=-00 

where is the pulse period. If the product x[n] = a[n]p[n] is filtered using a 
pulse form with impulse-response q[n\, then the result is the pulse-amplitude 
modulation(PAM) signal 

oO 

y[n]= Y^a\mN^\q\n-mN„\. (43) 

m=“Oo 

The cyclic spectra for the PAM signal, when a[n] is stationary, is given by 


s;(/)= 


0 , 


/ + X 
2J 


k 

To 

k 

T 


(44) 


where ^f) is the Fourier transform of the pulse form q[n], and S^f) is the 
power spectral density of the signal a[n]. 

Figures 24-29 show surface and contour plots of the SCD estimate 
magnitude for a PAM signal which is modulated by a random sequence of zeros 
and ones. 

From Eq. 44, ones expect to obtain something different from zero just at 
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cyclic frequencies that are multiple of the bit rate. Figures 24-29 verify the result 
obtained in Eq. 44. The poor resolution is more obvious in these examples, 
making it difficult to obtain the appropriate SCD representation for the signal. 
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Figure 24 Surface plot of the SCD estimate magnitude for a Pulse-Amplitude 
Modulated (PAM) Signal, using AUTOFAM, with the following parameters: 
A/ = 512/fe, Aa = 16/fe, rf^^XOlAHz, and /, =8192/fe, where is the bit rate, 
and /j is the sampling frequency. 
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Figure 25 Contour plot of the SCO estimate magnitude for a Pulse-Amplitude 
Modulated (PAM) Signal, using AUTOFAM, with the following parameters: 
A/ = 512ife, Aa = 16//z, r^=\02AHz, and /,=8192i/z, where is the bit rate, 
and is the sampling frequency. 
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Figure 26 Plots of the SCD estimate magnitude for a Pulse-Amplitude 
Modulated (PAM) Signal, using AUTOFAM, for a = 0, a = \Q2AHz, a = 204%Hz, 
a = 3012Hz, and a 4096Hz, respectively, with the following parameters: 
A/ = 512//z, Aa = 16i/z, /; = 1024/fe, and f, = m2Hz, where is the bit rate, 
and is the sampling frequency. 
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Figure 27 Surface plot of the SCD estimate magnitude for a Pulse-Amplitude 
Modulated (PAM) Signal, using AUTOSSCA, with the following parameters: 
A/ = 5\2Hz, Aa = \6Hz, =\024Hz, and =81927/2, where is the bit rate, 

and is the sampling frequency. 
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Figure 28 Contour plot of the SCD estimate magnitude for a Pulse-Amplitude 
Modulated (PAM) Signal, using AUTOSSCA, with the following parameters: 
A/" = 5\2Hz, A.a = \6Hz, =1024//z, and /^ = 8192/fe, where is the bit rate, 

and is the sampling frequency. 
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Figure 29 Plots of the SCD estimate magnitude for a Pulse-Amplitude 
Modulated (PAM) Signal, using AUTOSSCA, for a = 0 , or = 1024/fe , a = 2048//z, 
a = 3072/fe, and a = 4096//z, respectively, \with the following parameters; 
A/" = 512//z, Aa = 16/fe, r^ = 1024//z, and f^=Z\92Hz, where is the bit rate, 
and is the sampling frequency. 
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B. DIGITAL MODULATED SIGNALS 

1. Amplitude Shift Keying (ASK) Signal 

An ASK signal is simply an AM signal 


x\n\ = a\ri\ cos 


/ ^ 
Js 


(45) 


in which the amplitude message a[n] is a il/-ary PAM signal 


00 

a[n\ = - kN„ - . 


(46) 


A:=-oo 


The spectral correlation density function for this ASK signal is given by 






+e[/-/..+fj Q\f-f. -fj«:(/-/,) 


.-ilzarig 


+e(/+f+/,) -/,) 

+e[/+f-/•) e'(/-f+/.) sr=/-(/)e-'N-v.K-».i|, (47) 

For a full-duty cycle rectangle pulse, q[n] is given by 

* 1 = 6 : 

Its Fourier transform is given by 


/2 

/ 2 ' 


(48) 
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( 49 ) 



Details on this result can be found in Reference 9. 

Figures 30-41 present the outputs from AUTOFAM and AUTOSSCA for 
some ASK signals. Figures 30-35 present the results for an ASK signal with bit 
rate = 2048//z , and Figures 36-41 present the results for the same signal with 
/;=1024//z. As /„=2048//z, from Eq. 47, ones expect to obtain a 
representation which is a combination of four PAM signals centered at 
/ = = ±2048//z , for a = 0 , and at / = 0 , for a = ±2/, = ±4096//z . The results 


obtained in Figures 30-41 confirm that expectation. 

2. Binary-Phase Shift Keying (BPSK) Signal 

A phase-shift keying (PSK) signal is just a phase-modulated (PM) carrier 


f fo 

x[n] = cos 2n-^n-\-(j\n\ 
V Js 


(50) 


in which the phase (pin] is a M-ary PAM signal 


4«]= 'Y^a[mN^\q\n-mN„\. (51) 

m--<G 

Hence, the cyclic spectra for the BPSK signal, x[«], is given by Eq. 47, 
already computed for a M-ary ASK signal. Of course, the results are the same 
as for the M-ary ASK signal and are given by Figures 30-41. 
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Figure 30 Surface plot of the SCD estimate magnitude for an Amplitude Shift 
Keying (ASK) or a Binary Phase Shift Keying (BPSK) Signal, using AUTOFAM, 
with the following parameters: A/ = 512//z, Aa = 16//z, f^=204SHz, 

r^ = 204SHz, and f^ = S\92Hz, where is the carrier frequency, is the bit 
rate, and is the sampling frequency. 
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Figure 31 Contour plot of the SCO estimate magnitude for an Amplitude Shift 
Keying (ASK) or a Binary Phase Shift Keying (BPSK) Signal, using AUTOFAM, 
with the following parameters; A/ = 512/fe, Aa = 16^z, /^=2048i/z, 

r^=2QAZHz, and /^ = 8192^z, where is the carrier frequency, is the bit 
rate, and is the sampling frequency. 
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Figure 32 Plots of the SCD estimate magnitude for an Amplitude Shift Keying 
(ASK) or a Binary Phase Shift Keying (BPSK) Signal, using AUTOFAM, for 
a = 0, a = 2048//Z , a = 4096Hz , and a = 6144ffe , respectively, with the following 
parameters: A/ = 512^fe, La = \6Hz, f^=204%Hz, rj=2048//z, and 

/^ = 8192//z, where is the carrier frequency, is the bit rate, and is the 
sampling frequency. 
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Figure 33 Surface plot of the SCD estimate magnitude for an Amplitude Shift 
Keying (ASK) or a Binary Phase Shift Keying (BPSK) Signal, using AUTOSSCA, 
with the following parameters: A/ = 1024i/z, Aa = 8//z, f^=10A%Hz, 

r^=20A^Hz, and /^ = 8192//z, where is the carrier frequency, is the bit 
rate, and is the sampling frequency. 
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Figure 34 Contour plot of the SCD estimate magnitude for an Amplitude Shift 
Keying (ASK) or a Binary Phase Shift Keying (BPSK) Signal, using AUTOSSCA, 
with the following parameters: A/ = 1024/yz, Aa = 8ife, /^=2048//z, 
rj=2048//z, and /^ = 8192//z, where is the carrier frequency, is the bit 
rate, and is the sampling frequency. 
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Figure 35 Plots of the SCD estimate magnitude for an Amplitude Shift Keying 
(ASK) or a Binary Phase Shift Keying (BPSK) Signal, using AUTOSSCA, for 
a=0, a = 2048Hz, a = 4096Hz .and a = 6\44Hz , respectively, with the following 
parameters: A/ = 1024i/z, Aa = 8Hz, =2048Hz, r^= 2048Hz, and 

f^ = 8\92Hz, where is the carrier frequency, is the bit rate, and is the 
sampling frequency. 
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Figure 36 Surface plot of the SCD estimate magnitude for an Amplitude Shift 
Keying (ASK) or a Binary Phase Shift Keying (BPSK) Signal, using AUTOFAM, 
with the following parameters: A/ = 512//z, La = \6Hz, /.=2048^fe, 

r^=\(yiAHz, and f^ = Z\92Hz, where is the carrier frequency, is the bit 
rate, and is the sampling frequency. 
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Figure 37 Contour plot of the SCO estimate magnitude for an Amplitude Shift 
Keying (ASK) or a Binary Phase Shift Keying (BPSK) Signal, using AUTOFAM, 
with the following parameters: A/= 5127/2, Aa = 167/z, /, = 2048//z, 

r^ = \Q2AHz, and /^= 81927/2, where is the carrier frequency, is the bit 
rate, and is the sampling frequency. 
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Figure 38 Plots of the SCD estimate magnitude for an Amplitude Shift Keying 
(ASK) or a Binary Phase Shift Keying (BPSK) Signal, using AUTOFAM, for 
a = 0, a = 1024//z, a = 204SHz, a = 3072//z, and a = 4096//z, respectively, with 
the following parameters: A/ = 5\2Hz , Aa = \6Hz , /. = 2048/fz, = \024Hz , and 

f^ = S\92Hz, where is the carrier frequency, is the bit rate, and is the 
sampling frequency. 
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Figure 39 Surface plot of the SCD estimate magnitude for an Amplitude Shift 
Keying (ASK) or a Binary Phase Shift Keying (BPSK) Signal, using AUTOSSCA, 
with the following parameters; A/ = 512/fe, Aa = \6Hz, /^=2048/fe, 

?i=1024/fe, and f^ = S\92Hz, where is the carrier frequency, r* is the bit 
rate, and is the sampling frequency. 
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Figure 40 Contour plot of the SCD estimate magnitude for an Amplitude Shift 
Keying (ASK) or a Binary Phase Shift Keying (BPSK) Signal, using AUTOSSCA, 
with the following parameters; A/ = 5\2Hz , Aa = A Hz, /. = 2048^z , = 1024/fe, 

and /^ = 8192/fe, where is the carrier frequency, is the bit rate, and is 
the sampling frequency. 
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Figure 41 Plots of the SCD estimate magnitude for an Amplitude Shift Keying 
(ASK) or a Binary Phase Shift Keying (BPSK) Signal, using AUTOSSCA, for 
a = 0, a = 1024Hz, a = 2048/fe, a = 3072/fe , and a = 4096Hz , respectively, with 
the following parameters: A/ = 5l2Hz , Aa = 16Hz , /<, = 2048/fe , = \024Hz , and 

= 8192/fe, where is the carrier frequency, is the bit rate, and f, is the 
sampling frequency. 
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VI. CONCLUSIONS 


A. SUMMARY 

The purpose of this thesis was to implement the FAM and the SSCA 
methods in MATLAB, allowing students, researchers, and engineers to take 
advantage of the power of the cyclic analysis methods for solving signal 
detection and modulation identification problems. Since MATLAB is user friendly 
and easily portable to other operating systems, the implementation becomes a 
proving ground, easily modified and set up on other computer systems. 

The two methods are implemented as user friendly as possible. The only 
inputs required for both are the signal(s), the sampling frequency (/J, that 
should be the same for both, and the desired resolutions for spectrum frequency 
(/) and cyclic frequency (a). The results generated by both programs can be 
displayed in different ways, such as surface plots, contour plots, and cross- 
section plots. 

Both programs generate a large amount of data when a good resolution is 
desired. As a consequence, limitations in computational as well as printing 
resources did not allow the presentation of more detailed plots. One can of 
course easily create averaged spectral correlation surfaces, with a coherent or 
power averaging method. 
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B. SUGGESTIONS 


A step that was not possible to reach during this work was an analysis 
and experimentation of the performance of cyclic spectral analysis in a white- 
noise contaminated environment. A natural extension to this work is an analysis 
of how well cyclic analysis performs in a white and a colored noise environment 
such as the one imposed by the ocean to the receiving elements. 

Computationally, cyclic spectral analysis is an expensive task. Therefore, 
any improvements to the speed of existing methods or even completely new fast 
algorithms are desirable. 

During the development of this thesis, it was intended to replace some of 
the Fast Fourier Transform (FFT) operations by the fast Wavelet Transform. A 
follow on to this work is the replacement of the FFT operations by Wavelet 
Transforms, and an evaluation of how well the modified algorithm performs in 
terms of identification and computational cost. 
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APPENDIX A. CALCULATION OF THE SCD FUNCTION OF AN 
AMPLITUDE-MODULATED SIGNAL 


Let us consider the amplitude-modulated (AM) signal 

x[«] = a[n\ cos(2;r / , (52) 

where a[n] is a stationary random lowpass signal with PSD 5„(/), with no 
spectral lines in its PSD. 

The fundamental parameter for second-order periodicity in discrete-time, 
called cyclic autocorrelation function and denoted by . is given by 

[£] = {x[«] x'[n- ) e""'. (53) 


We can look at R^[£] from the following point of view; 

Let u[n\ = x[«] and v[«] = x[n\ , then 

K ¥\={^[«] e'"" [« - 
= [« - ) 

= («[«] v*[«- ^]) 

= ^.[4 ( 54 ) 

But, since u[n\ = x[«] , 
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= (“[«] u[n-i^ 

= {x[«]x'[«-4e''^ 

= (55) 

and taking the Fourier transform leads to 

5.(/) = S,(/)®5(/ + f]=S.(/ + |). (56) 

Similarly, if v[«] = x[n] e““", then 

r , r , 

«.W=*,W« ' . (57) 

and 

S.(/) = S,(/-f). (58) 

Thus, we can redefine the second-order periodicity, by saying that x[n\ 
exhibits second-order periodicity if and only if frequency translates (frequency- 
shifted versions) of x[n] are correlated with each other. 

As long as the mean values of the frequency translates u[n\ and v[n\ are 
zero (i. e., («[«]) = 0 and (v[«]) = 0), x[n\ does not contain finite-amplitude 
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additive sine-wave components at ±a/2 and therefore -S',(/) has no spectral 
lines at / = ±a/2. The crosscorrelation is actually a temporal 


crosscovariance K^[i], that is, 


= ({«[«] - (“[«])} {v[« - -f] - (v[« - f])} 
= («[«] v*[n-^]) 




(59) 


Based on this, we can define a temporal correlation coefficient, the 
magnitude of which is upper bounded by unity, for frequency translates since 


■''.loj (/:,[0]x.[0])’ 


(60) 


Since a[n] is a real stationary random ( R“[£] = 0, for all a^O) signal with 
zero mean (i. e., (a[n]) = 0 ), the cyclic autocorrelation function of x[«] is given 


by 


K W = (^[«] x[n-£] ) e'" 


i7ra£ 


= (a\n]cos(2rr f^n + O) a‘[n-£]cos[2;r/„(«-£)+^]e 

= /a[n] a* \n - £]— |cos(2;7’ /,,£) + cos[2;r/„ (2n -£)+ 2^]| \ 
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= I cos(2;r//) {a[n\ a[n- 

+ ^(a[/i]a* [« - f]cos[2;r /o(2« - £)+ 

1 1 / ^42^/o( 2«-^)+2^] -/[2;r/^(2n-£)+2^] 


+ ic"“(o[H]a' [n - <] ) jM**;/.)' 

=^^:WM2''//)+y'’" C^'W+y"'" 

( 61 ) 

But since a[n] is a stationary signal (i. e., R^[l] = 0 for all a * 0), the only 
non-zero contributions are at a = 0 and a = ±2f„. Thus, the cyclic 
autocorrelation function becomes 


KU] = 


|/?J^]cos(2;r//), 

0 , 


a = ^fo 

« = 0 
otherwise 


( 62 ) 


Defining the temporal correlation coefficient as 
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(63) 


s:[£]= 




where R^[0] = R“[£] for a = 0, computed at ^ = 0, so R^[o] = i?„[0]/2 
So, S°[l] becomes 


s:[e]= 






= —e 


±ne 


m 


(x = ±2fo 


- = f,t), a = 0. 




(64) 


Thus, the strength of correlation between and x[«-f] is 


given by 


k‘WI= 


a = ±2f, 

|<5„[<]cos(2;>r/,£)|, a = 0 
0, otherwise. 


(65) 


It can be substantial for an amplitude-modulated signal, e.g.,|<?“[0]| = 1/2. 

To localize in the frequency domain the average power (|x[n]p) = [O] in 

a stationary signal x[n\, we simply pass the signal x[«] through a set of 
narrowband bandpass filters and then measure the average power at the output 
of the filters. In the limit when the bandwidths of these filters approach zero, the 
corresponding set of measurements of average power, normalized by the 
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bandwidth, constitute the "power spectral density ( PSD ) function", given by 

‘5.(/) = f), (66) 

where hl[n\ is the discrete-impulse response of a bandpass filter with center 
frequency /, bandwidth B , and unity gain at the band center. 

In a similar fashion, to localize in the frequency domain the correlation 

(m[«]v’[«]) = =/?“[0] of frequency-shifted versions u[n] and v[«] of 

a cyclostationary signal x[n], we pass both of the two frequency translates u[n] 
and v[n] of x[«] through the same set of bandpass filters. Again, in the limit 
when the bandwidth of these filters approach zero, the corresponding set of 
measurements of the temporal correlation of the filtered signals constitute the so 
called "spectral-correlation density (SCD ) function", given by 

(/) = W[^1 ® [«] ® v[«]} ^, (67) 

that yields the spectral density of correlation between u[n] and v[n] at frequency 
/, which is identical to the spectral density of correlation in x[n] at frequencies 
/ + a/2 and / - a/2. 

It's a well known result, called "Wiener relation" as opposed to its 
probabilistic counterpart called "Wiener-Khinchin relation", that the PSD is equal 
to the Fourier transform of the autocorrelation function. 
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( 68 ) 


«,(/)= I 

^=-00 

Similarly, the SCD can be obtained by Fourier transforming the "cyclic 
autocorrelation function", 


^.'(/)= Z (69) 

£~—co 

This result is known as the "cyclic ( periodic) Wiener relation". 

We observe that since R“[£] is periodic in a with period two, so is 

s;{f) 

[^] = x'{n-£] ^ 

= (x[n] X* [n - i] , (70) 

and, since n and £ are integers, 

= 1, (71) 


so. 


¥] = W x'{n-£] 

= K[£]- 


(72) 


Thus, 
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sr(/)= E 


-ilrr Ji 


^=-00 


= Z 

= «,“(/). (73) 


Also, since i takes on only integer values, then S“{f) is periodic in / 
with perio(d one 


5; (/+1) = Z K [^] 


^=—00 


= E 


i=-<o 


e-n^ 


(74) 


and, since £ is an integer, we have 


^-,2^ 


= 1 . 


(75) 


then. 


s; (/+!)= E 

i=-co 

Furthermore, 5“(/) also exhibits the following periodicity: 


(76) 


sr'(/ + i/2) = s;(/). 


(77) 
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This is easily shown by 




^=-00 




-iTti 


i=-oo 


= 'Z{x[n]x[n-e]e-‘^’^‘^ e 


-iTii 


£=-00 




£=-00 


(78) 


and, since n is an integer, we have 


e-'^’^ = 1, 


(79) 


hence. 


5 ;" [/ += X{x[rt]x* [n - £] e"'""”) 6““" 


£=—00 




£=-00 


=5r(/)- 


(80) 


We also define a spectral correlation coefficient, p“{f ), given by 


A“(/) = 


«,“(/) 
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(81) 


S..(/) 


[s,(/K(/)]^ 


Since is bounded by [0,l], it is a convenient measure of the 

degree of local spectral redundancy that results from spectral correlation. 

Going back to the AM signal, by Fourier transforming the cyclic 
autocorrelation function, we get 


sr(/)= 


a=±2/. 


0 , 


otherwise 


(82) 


and so, the spectral correlation coefficient is given by: 


/ + yJ = 5;(/), for a = 0, computed at Z + y 

= + fo) + \sXf - fo)^ computed at / + ^ 




(83) 


Also, 




(84) 


So that, fora = 2Z,, 
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■s, (/ + fJ = s, (/ + /,) = i S, (/ + 2/„) + i S, (/), ( 85 ) 

and, 

■s- [/ - fJ = s. (/ - /. ) = J (/) + J 5, (/ - 2/.). (86) 

and for « = -2/^, 

s, (/ + f) = S, (/ - /.) = J s, (/) + J s.(/ - 2/,), (87) 

and, 

5- [/ - fJ = s. (/ + /,) = J s, (/ + 2/,) + J S, (/), (88) 

and, for a = 0, 

s, (/ + f) = s, (/) = j s, (/ + /.) + ^5, (/-/,), (89) 

and, 

■S',(/ - f) = S,(/) = jS.(/ + /.) + ^s.(/- /.), (90) 



sAfV"" 

II 

1 

{K (/ + 2/J + 5, (/)][5,(/) + sXf- )]p 

1, 

0, 


« = ± 2 /„ 


a = 0 
otherwise. 


(91) 
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Thus, the strength of correlation between the spectral components of x{n] 
at frequencies / + a/2 and / -a/2 is unity (i. e., |/?“(/)| = l, for |/| </„ and 
a = ±2/„), provided that a\n] is bandlimited to l/|</<,. Hence 5„(/) = 0 for 

I/I a f. ■ 
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APPENDIX B. FUNCTION AUTOFAM 


function [Sx, alphao, fo] =autofatn (x, f s, df, dalpha) 

%■ AUTOFAM (X, FS, DF, DALPHA) computes the spectral auto-correlation 

% density function estimate of the signal X, by using the FFT 

% Accumulation Method (FAM) . Make sure that DF is much bigger 

% than DALPHA in order to have a reliable estimate. 


INPUTS: 

X 

FS 

DF 

DALPHA - 


input column vector; 
sampling rate; 

desired frequency resolution; and 
desired cyclic frequency resolution. 


OUTPUTS: 

SX - spectral correlation density function estimate; 

ALPHAO - cyclic frequency; and 
FO - spectrum frequency. 


Author: E.L.Da Costa,9/28/95. 


if nargin ~= 4 

error('Wrong number of arguments.'); 

end 


%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% Definition of Parameters % 

%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

Np=pow2(nextpow2(fs/df)); % Number of input channels, defined 

% by the desired frequency 
% resolution(df) as follows: 

% Np=fs/df, where fs is the original 
% data sampling rate. It must be a 
% power of 2 to avoid truncation or 
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% zero-padding in the FFT routines; 

L=Np/4; % Offset between points in the same 

% column at consecutive rows in the 
% same channelization matrix. It 
% should be chosen to be less than 
% or equal to Np/4; 

P=pow2(nextpow2(fs/dalpha/L)); % Number of rows formed in the 

% channelization matrix, defined 
% by the desired cyclic frequency 
% resolution(dalpha) as follows: 

% P=fs/dalpha/L.It must be a power 
% of 2/ 

N=P*L; % Total number of points in the 

% input data. 




Input Channelization 


9-9-9-9-9-9-9-9-9-2.9-9-S'S'9'S'5'9-S-e'S-9-S'9 
oo ooo 0^0 o^^^o'S’ooo^^ t 


if length(x)<N 
x(N)=0; 

elseif length(x)>N 
x=x(1:N); 

end 


NN= (P-1) *L+Np; 


xx=x; 
xx(NN)=0; 
xx=xx(:); 

X=zeros(Np,P); 
for k=0:P-1 

X(:,k+l)=xx(k*L+l:k*L+Np); 


end 


S-S-5-9-9-9-9-S-9-5>2-9-S- 

% Windowing % 
%%%%%%%%%%%%% 
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a=hamming (Np) ; 
XW=diag{a)*X; 
XW=X; 


O'O'O'O'O'O O 0*0 o'o'o o 

% First FFT % 

2'9-2'S-S'2'9-2'2'S-9'9-9- 

ooooooooooooo 

XFl=fft(XW); 

XFl=fftshift(XFl); 

XF1=[XF1{:,P/2+l:P) XFl{:,1:P/2)]; 


S-9.9.2-S.S.9.S.2-2-9-e-2-9'9-9-e-2- 

000000000000 * 00000*0 

% Downconversion % 


E=zeros(Np,P); 
for k=-Np/2:Np/2-l 
for m=0:P“l 

E(k+Np/2+1,m+1)=exp(-i*2*pi*k*m*L/Np) ; 

end 

end 


XD=XF1.*E; 
XD=conj(XD'); 


% Multiplication % 

S'S'9-2-9-9-S-9-9-2-9-2-2-9-S-S-9-9- 

o^*o*o*o*o*o*o*0'o*o*o‘o*o*o*o*o*o 

XM=zeros(P,Np^2); 
for k=l:Np 

for 1=1:Np 

XM{:,(k-1)*Np+l)={XD{:,k). *COnj (XD{:,1))) 

end 

end 
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% Second FFT % 

9-9-S-S-2-9-S-2-2-2-2-2-2-2- 
'o'o'o'o'o'o'o'o'o o'o'o o o 

XF2=fft(XM); 

XF2=fftshift(XF2); 

XF2=[XF2(:,Np^2/2+l:Np^2) XF2(:,l:Np^2/2)] 
XF2=XF2(P/4:3*P/4,:); 

M=abs{XF2); 
alphao=-l:1/N:1; 
fo=-.5:l/Np:.5; 

Sx=zeros(Np+1,2*N+1); 
for kl=l:P/2+l 

for k2=l:Np^2 

if rem(k2,Np)==0 
l=Np/2-l; 
else 

l=rem(k2,Np)-Np/2-1; 

end 

k=ceil(k2/Np)-Np/2-1; 
p=kl-P/4-l; 

alpha=(k-1)/Np+(p-1)/L/P; 
f=(k+l)/2/Np; 
if alpha<-l | alpha>l 
k2=k2+l; 

elseif f<-.5 | f>.5 
k2=k2+l; 
else 

kk=l+Np*(f+.5); 

11=1+N*(alpha+1); 

Sx(kk,ll)=M(kl,k2); 

end 

end 

end 
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APPENDIX C. FUNCTION AUTOSSCA 


function [Sx,alphao,fo]=autossca(x,fs,df,dalpha) 

% AUTOSSCA(X,FS,DF,DALPHA) computes the spectral auto- 

% correlation density function estimate of the signal X, 

% by using the Strip Spectral Correlation Algorithm (SSCA). 

% Make sure that DF is much bigger than DALPHA in order to 

% have a reliable estimate. 


INPUTS: 

X 

FS 

DF 

DALPHA - 


input column vector; 
sampling rate; 

desired frequency resolution; and 
desired cyclic frequency resolution. 


OUTPUTS: 

SX - spectral auto-correlation density function estimate; 

ALPHAO - cyclic frequency; and 
FO - spectrum frequency. 


Author: E.L.Da Costa,9/28/95. 


if nargin -= 4 

error('Wrong number of arguments'); 

end 


% Definition of Parameters % 


oooooooooooooooooooooooooooo 


Np=pow2(nextpow2(fs/df)); 


% Number of input channels, defined 
% by the desired frequency 
% resolution(df) as follows: 

% Np=fs/df, where fs is the original 
% data sampling rate. It must be a 
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L=Np/4; 


P=pow2(nextpow2(fs/dalpha/L))/ 


N=P*L; 


% power of 2 to avoid trxincation or 
% zero-padding in the FFT routines; 
% Offset between points in the same 
% column at consecutive rows in the 
% same channelization matrix. It 
% should be chosen to be less than 
% or equal to Np/4; 

% Number of rows formed in the 
% channelization matrix, defined by 
% the desired cyclic frequency 
% resolution(dalpha) as follows: 

% P=fs/dalpha/L. It must be a power 
% of 2; 

% Total number of points in the 
% input data. 


% Input Channelization % 


if length(x)<N 
x(N)=0; 

disp('you will not get the desired resolution in cyclic frequency'); 
dalpha=fs/N; 

disp(['cyclic frequency resolution=',num2str(dalpha)]); 
elseif length(x)>N 
x=x(1:N); 


end 

NN= (P-l) *L+Np; 


XX=X; 

XX(NN)=0; 
xx=xx(:); 

X=zeros(Np,P); 
for k=0:P-l 

X ( :,k+1)=xx{k*L+l:k*L+Np) / 

end 
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% Windowing % 


a=hamming(Np); 
XW=diag(a)*X; 


S'S'2-e-S-2-S-2'2-2-S-2-2- 

% First FFT % 


i-e-s-s-s-s-s-s-s-g- 


XFl=fft(XW); 

XFl=fftshift(XFl); 

XF1=[XFl(;,P/2+1:P) XFl(:,1:P/2)]; 


g'2-g-S-g-2-&-2-&-8-e-2-B-S-e-9-S-e- 
o o o o'o'p'o'o'o'o'o'o'o'o'o'o'o'b 

% Downconversion % 


>00000000000000000 


E= 2 eros(Np,P); 
for k=-Np/2:Np/2-l 
for m=0:P-1 

E(k+Np/2+l,m+l)=exp(-i*2*pi*k*m*L/Np) 

end 


end 


XD=XF1.*E; 


000000000000000 

% Replication % 

S-9-S-9-S-2-S-9-a-e-9-S.Q-S.9- 

■o'o'o'o'b'o'o'o'o^'o'o'o'o'o 

XR=2eros(Np,P*L); 
for k=l;P 

XR(:,(k-l)*L+l:k*L)=XD(:,k)*ones(l,L); 

end 


% Multiplication % 
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oooooooooooooooooo 


xc=ones(Np,1)*x'; 
XM=XR.*xc; 

XM=conj(XM'); 


% Second FFT % 


XF2=fft(XM); 

XF2=fftshift(XF2); 

XF2=[XF2(:,Np/2+l;Np) XF2(:,1:Np/2)]; 

M=abs(XF2); 

alphao=(-1:l/N:1)*fs; 

fo=(-.5:1/Np:.5)*fs; 

Sx=zeros(Np+1,2*N+1); 
for kl=l:N 


for k2=l:Np 

alpha=(kl-1)/N+(k2-l)/Np-1; 
f={(k2-l)/Np-(kl-1)/N)/2; 
k=l+Np*(f+.5); 

1=1+N*(alpha+1); 
Sx(k,l)=M(kl,k2); 

end 


end 
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APPENDIX D. FUNCTION CROSSFAM 


function [Sxy,alphao,fo]=crossfam(x,y,fs,df,dalpha) 

% CROSSFAM(X,Y,FS,DF,DALPHA) computes the spectral cross- 

% correlation density function estimate of the signals X 

% and Y, by using the FFT Accumulation Method(FAM). Make 

% sure that DF is much bigger than DALPHA in order to have 

% a reliable estimate. 


INPUTS : 

X 

Y 

FS 

DF 

DALPHA - 


input column vector; 
input column vector; 
sampling rate; 

desired frequency resolution; and 
desired cyclic frequency resolution. 


OUTPUTS: 

SXY - spectral cross-correlation density function estimate; 
ALPHAO - cyclic frequency; and 
FO - spectrum frequency. 


Author: E.L.Da Costa,9/28/95. 


If nargin ~= 5 

error('Wrong number of arguments.'); 

end 


% Definition of Parameters 


oooooooooooooooo 


Np=pow2(nextpow2(fs/df)); 


% Number of input channels, defined 
% by the desired frequency 
% resolution(df) as follows: 

% Np=fs/df, where fs is the original 
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L=Np/4; 


P=pow2(nextpow2(fs/dalpha/L)); 


N=P*L/ 


% data sampling rate. It must be a 
% power of 2 to avoid truncation or 
% zero-padding in the FFT routines; 
% Offset between points in the same 
% column at consecutive rows in the 
% same channelization matrix. It 
% should be chosen to be less than 
% or equal to Np/4; 

% Number of rows formed in the 
% channelization matrix, defined by 
% the desired cyclic frequency 
% resolution(dalpha) as follows: 

% P=fs/dalpha/L. It must be a power 
% of 2 ; 

% Total number of points in the 
% input data. 




% Input Channelization % 

if length(x)<N 
x(N)=0; 

elseif length(x)>N 
x=x (1 :N) ; 

end 

if length(y)<N 
y(N)=0; 

elseif length(y)>N 
y=y(1:N); 

end 


NN= (P-1) *L+Np; 


xx=x; 


yy=y; 

XX (NN)=0; 
yy(NN)=0; 
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xx=xx{:); 
yy=yy(;); 

X=zeros{Np,P)/ 

Y=zeros{Np,P); 
for k=0:P-1 

X ( : , k+1) =xx (k*L+l :k*L+Np) ; 
Y ( : , k+1) =yy (k*L+l :k*L+Np) ; 

end 


S-S-S-S-S-S-S-S-S-S-S-S-S- 

oo'^oooooooooo 

% Windowing % 


a=hamming (Np) / 
XW=diag(a)*X; 
YW=diag(a)*Y; 


2'S'2-S'S-S-2-9'9-2-S-S-9- 

0000000*000000 

% First FFT % 

%%%%%%%%%%%%% 

XFl=fft(XW); 

YFl=fft(YW); 

XFl=fftshift(XFl); 

YFl=fftshift(YFl); 

XF1=[XFl{:,P/2 + 1:P) XFl{:,l:P/2) ] ; 
YF1=[YFl(:,P/2 + 1:P) YFl{:,l:P/2) ] ; 


oooooooooooooooooo 


% Downconversion % 


S'S-S'S-S-S-S'S'S'S'S'S'S'S'S'S'S'S- 


E=zeros(Np,P); 
for k=-Np/2+l:Np/2 
for m=0:P-l 

E {k+Np/2 ,m+l) =exp {-i*2*pi*k*m*L/Np) 


end 


end 
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XD=XF1.*E; 
YD=YF1.*E; 
XD=conj(XD'); 
YD=conj(YD'); 


2'2'9-8'9'9-9'2'9-S-9'S^S'9'9'9-9-2' 

00*6 o o'o'oo'oo'o o'o'o'b o o o 

% Multiplication % 

%%%%%%%%%%%%■%%%%%% 

XYM=zeros(P,Np^2); 
for k=l:Np 

for 1=1:Np 

XYM(:,(k-1)*Np+l)=XD(:,k).*conj(YD{:,1)) 

end 

end 

%%%%%%%%%%%%%% 

% Second FFT % 

2'2'S'5-9-2'9-9-S'S«S'9'2'9- 

oooooooooooooo 

XYF2=fft(XYM); 

XYF2=fftshift(XYF2); 

XYF2=[XYF2(:,Np^2/2+l:Np^2) XYF2(:,l:Np^2/2)]; 
XYF2=XYF2(P/4:3*P/4,:); 

M=abs(XYF2); 
alphao=(-1:l/N:1)*fs; 
fo=(-.5:1/Np:.5)*fs; 

Sxy= 2 eros(Np+1,2*N+1); 
for kl=l:P/2+l 

for k2=l:Np^2 

if rem(k2,Np)==0 
l=Np/2; 
else 

l=rem(k2,Np)-Np/2; 

end 

k=ceil(k2/Np)-Np/2; 
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p=kl-P/4-l; 

alpha=(k-1)/Np+(p-1)/L/P; 
f=(k+l)/2/Np; 
if alpha<>l | alpha>l 
k2=k2+l; 

elseif f<-.5 | f>.5 
k2=k2+l; 
else 

kk=l+Np*(f+.5); 

11=1+N*(alpha+1); 

Sxy(kk,11)=M(kl,k2); 

end 

end 

end 
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APPENDIX E. FUNCTION CROSSSSCA 


function [Sxy,alphao,fo]=crossssca(x,y,fs,df,dalpha) 

% CROSSSSCA(X,y,FS,DF,DALPHA) computes the spectral cross- 

% correlation density function estimate of the signals X 

% and Y, by using the Strip Spectral Correlation Algorithm 

% (SSCA) . Make sure that DF is much bigger than DALPHA in order 

% to have a reliable estimate. 

% 


% 

% 


INPUTS: 

X 

Y 

FS 

DF 

DALPHA - 


input column vector; 
input column vector; 
sampling rate; 

desired frequency resolution; and 
desired cyclic frequency resolution. 


OUTPUTS: 

SXY - spectral cross-correlation density function estimate; 
ALPHAO - cyclic frequency; and 
FO - spectrum frequency. 


Author: E.L.Da Costa,9/28/95. 


If nargin 5 

error('Wrong number of arguments.'); 

end 


S-9-S.S-2>9-S'9-9'9-S'2-9-2-9'9-9'9-9.9-2>9-9'9-9-9-9-9- 

oo'^'S'Sooooooooooooo'S'S'Soooooo'S 


% Definition of Parameters % 


9-9'2-9-9-9-&'S'9-9-9-9'9-9-9-9-9-9-9'9-9-9-9-9-9-9-9-9- 

oooooooooooooooooooooooooooo 


Np=pow2(nextpow2(fs/df) ) ; 


% Number of input channels, defined 
% by the desired frequency 
% resolution(df) as follows: 

% Np=fs/df, where fs is the original 
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% data sampling rate. It must be a 
% power of 2 to avoid truncation or 
% zero-padding in the FFT routines; 

L=Np/4; % Offset between points in the same 

% column at consecutive rows in the 
% same channelization matrix. It 
% should be chosen to be less than 
% or equal to Np/4; 

P=pow2(nextpow2(fs/dalpha/L)); % Number of rows formed in the 

% channelization matrix, defined by 
% the desired cyclic frequency 
% resolution(dalpha) as follows: 

% P=fs/dalpha/L. It must be a power 
% of 2; 

N=P*L; % Total number of points in the 

% input data. 


% Input Channelization % 

oooooooooooooooooooooooo 

if length{x)<N 
x{N)=0; 

elseif length(x)>N 
x=x (1 :N) / 

end 

if length(y)<N 
y{N)=0; 

elseif length(y)>N 
y=y {1:N) ; 

end 

NN=(P-1)*L+Np; 

XX=X; 

XX (NN) =0 ; 

xx=xx(:)/ 

X=zeros(Np,P); 
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for k=0:P-1 

X {:, k+1) =XX {k*L+l:k*L+Np); 

end 


% S-9'9-9'9'9-e-9'9'S-a-S- 
'o'o o o o 

% Windowing % 

S'9'9'5-5'S'2-S'S'S'S-2'&' 

00000000*00000 

a=hamming (Np) ; 

XW=diag(a)*X; 


% First FFT % 


2r9'2'2»2-9'St»2'2-2'St'2'2' 

"o o'o'o'o o'o'o'o'o'o'o'o 

XFl=fft(XW); 

XFl=fftshift(XFl); 

XF1=[XF1(:,P/2+l:P) XFl(:,1:P/2)]; 


%%%%%%%%%%%%%%%%%% 

% Downconversion % 

%%%%%%%%%%%%%%%%%% 

E=zeros(Np,P); 
for k=-Np/2+l:Np/2 
for m=0:P-l 

E(k+Np/2,m+l)=exp(-i*2*pi*k*m*L/Np) 

end 

end 

XD=XF1.*E; 

XD=conj(XD*); 

‘o'o'o'O'O^'O'o'o'o'o'S'o'o'o 

% Replication % 

9-9'2-S'a-e-S'S-9.S-9.S-e-9-S> 


XR=zeros(Np,P*L); 
for k=l;P 
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XR(:,(k-1)*L+1:k*L)=XD(:,k)*ones(1,L); 


end 


% Multiplication % 

'o'o’o o 0'O‘O'O'O‘S‘O‘O‘S o'o'o’o'o 

yc=ones(Np,1)*y'; 
XYM=XR. *yc; 
XYM=COnj(XYM'); 


% Second FFT 


8-S-S-9-2-9-9-9-9-9-a-a-S-9- 
o o o'o'o'o’o'o'o'o'o'o 0*0 


Second FFT 


XYF2=fft(XYM); 

XYF2=fftshift(XYF2); 

XYF2=[XYF2(:,Np"2/2+l:Np"2) XYF2(:, 
M=abs(XYF2); 
alphao=(-l:l/N:l)*fs; 
fo=(-.5:l/Np;.5)*fs; 

Sxy=zeros(Np+1,2*N+1); 
for kl=l:N 

for k2=l:Np 

alpha=(kl-l)/N+(k2-l)/Np-1; 
f=((k2-l)/Np-(kl-l)/N)/2; 
k=l+Np*(f+.5); 

1=1+N* (alplia+1) ; 

Sxy (k, 1) =M (kl, k2) ; 

end 

end 


l:Np"2/2)]; 


I 
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APPENDIX F. PLOTTING ROUTINES 


%Surface Ploti 

g-9-9-9-S'S-5-2-9- 

figure(1) 

surf1(alphao,fo,Sx); 
view(-37.5,60); 

title('SCD estimate using FAM'); 
xlabel('alpha'); 
ylabel('f'); 
zlabel(’Sx’); 

%%%%%%%%%%%%%% 

%Contour Plot% 

%%%%%%%%%%%%%% 
figure(2) 

contour(alphao,fo,Sx); 
xlabel('alpha(Hz)'); 
ylabel('f(Hz)'); 


%Cross-Section Plots% 


figure(3) 

plot(fo,Sx(:,l+N*(alpha/fs+1)) 

xlabel('f(Hz)'); 
ylabel('Sx(alpha)'); 


% alpha is the desired cyclic 
% frequency. 
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